 functions {
  real pct_correct_pred(int[] y, vector mu) {
    real out;
    int N;
    N = num_elements(mu);
    out = 0.;
    for (n in 1:N) {
      if (y[n]) {
        out = out + int_step(mu[n] >= 0.5);
      } else {
        out = out + int_step(mu[n] < 0.5);
      }
    }
    out = out / N;
    return out;
  }
  real expected_pct_correct_pred(int[] y, vector mu) {
    real out;
    int N;
    N = num_elements(mu);
    out = 0.;
    for (n in 1:N) {
      if (y[n]) {
        out = out + mu[n];
      } else {
        out = out + (1. - mu[n]);
      }
    }
    out = out / N;
    return out;
  }
}



data{
int<lower=1> N; //number of observations
int<lower=1> C; //number of oblasts
int<lower=1> K; //nubmer of covariates
int<lower=1> P; //number of parites
int<lower=1, upper=C> city[N];
int<lower=1, upper=P> party[N];
int<lower=0, upper=1> y[N];
matrix[N,K] x;
}


parameters{
real alpha;
vector[K] beta;
vector[C] ranef_c;
vector[P] ranef_p;
real<lower=0, upper=100> sigma_c;
real<lower=0, upper=100> sigma_p;
real<lower=3, upper=7> nu;
real<lower=3, upper=7>nu_c;
real<lower=3, upper=7>nu_p;
}

transformed parameters {
vector[N] theta;
vector[N] y_hat;
theta=x*beta;


for(n in 1:N){
y_hat[n]=alpha + theta[n] + ranef_p[party[n]] + ranef_c[city[n]];
}
}

model{
//priors
target+=normal_lpdf(alpha|0,1);

for(k in 1:K){
target+=student_t_lpdf(beta[k]| nu,0,2.5);
}
target+=student_t_lpdf(sigma_c|nu_c,0,1);
target+=student_t_lpdf(sigma_p|nu_p,0,1);
target+=gamma_lpdf(nu|2, 0.1);
target+=gamma_lpdf(nu_p|2, 0.1);
target+=gamma_lpdf(nu_c|2, 0.1);

to_vector(ranef_c)~normal(0,sigma_c);
to_vector(ranef_p)~normal(0,sigma_p);

//likelihood
target+= bernoulli_logit_lpmf(y|y_hat);



}

generated quantities{
int y_rep[N];
vector[N] log_lik;
real PCP;
real ePCP;
vector[N] mu;
vector[K] fd_beta;


for(n in 1:N){
  mu[n]=inv_logit(y_hat[n]);
  log_lik[n]=bernoulli_logit_lpmf(y[n]|y_hat[n]);
  y_rep[n]=bernoulli_logit_rng(y_hat[n]);
}

PCP=pct_correct_pred(y, mu);
ePCP=expected_pct_correct_pred(y,mu);


for(k in 1:K){
  if(k==4)
  fd_beta[k]=exp(alpha +beta[k]*1)/(1+exp(alpha +beta[k]*1))-exp(alpha +beta[k]*0)/(1+exp(alpha +beta[k]*0));
  else if(k==5)
  fd_beta[k]=exp(alpha +beta[k]*1)/(1+exp(alpha +beta[k]*1))-exp(alpha +beta[k]*0)/(1+exp(alpha +beta[k]*0));
  else
  fd_beta[k]=exp(alpha +beta[k]*(mean(x[,k]+sd(x[,k]))))/(1+exp(alpha +beta[k]*(mean(x[,k]+sd(x[,k])))))-exp(alpha +beta[k]*(mean(x[,k]-sd(x[,k]))))/(1+exp(alpha +beta[k]*(mean(x[,k]-sd(x[,k])))));
}


}
